Topological multicritical point in the Toric Code and 3D gauge Higgs Models 
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We report a new type of multicritical point that arises from competition between the Higgs and 
confinement transitions in a Z2 gauge system. The phase diagram of the 3d gauge Higgs model has 
been obtained by Monte-Carlo simulation on large (up to 60"^) lattices. We find the transition lines 
continue as 2nd-order until merging into a Ist-order line. These findings pose the question of an 
effective field theory for a multicritical point involving noncommuting order parameters. A similar 
phase diagram is predicted for the 2-dimensional quantum toric code model with two external fields, 
hz and hx\ this problem can be mapped onto an anisotropic 3D gauge Higgs model. 

PACS numbers: 



Introduction. Topological quantum phases and anyons 
are well known in connection with the fractional quantum 
Hall effect, but they are also expected to exist in frus- 
trated magnets. It has long been proposed that a certain 
class of resonating- valence-bond (RVB) [l[ phases carries 
Z2-charges and vortices and has a four-fold degenerate 
ground state on a torus Q]. A qualitative understand- 
ing of this phase can be obtained from a so-called toric 
code model (TCM) The dimer model on the Kagome 
lattice is mapped onto the TCM exactly Q while some 
other models [1,0] belong to the same universality class. 

The TCM is defined in terms of spin- 1/2 degrees of 
freedom located on the bonds of an arbitrary 2d lattice. 
The Hamiltonian is as follows: 
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and Bp = Hjep "^1 products of 
spin operators (u" arc the Pauli matrixes) on the bonds 
incident to a lattice site s and on the boundary of a pla- 
quette p, respectively. The ground state corresponds to 
eigenvalues Ag ~ \, Bp = 1 for all s and p. On a surface 
of genus g, it is 4^-fold degenerate. Elementary excita- 
tions are characterized by eigenvalues Ag = —1 (a Z2 
charge on site s) and Bp = — 1 (a Z2 vortex on plaquette 
p); all excitations are gapped. Each type of quasiparticle 
is bosonic, but due to nontrivial mutual braiding, they 
must be jointly regarded as Abelian anyons. 

The Hamiltonian ([T|) has special properties related to 
its exact solvability: the two-point correlator vanishes 
and the quasiparticles have flat dispersion. These fea- 
tures do not survive a small generic perturbation, while 
the topological character of the ground state and the any- 
onic quasiparticle statistics are robust. Yet a sufficiently 
strong field can polarize the spins, driving a transition 
to the topologically trivial phase. Trebst at al stud- 
ied a perturbation of the form —h and solved the 
problem by reducing it to the 2d transverse-field Ising 



model, which is mapped to an anisotropic 3d classical 
Ising model. In this paper we consider a more general 
Hamiltonian: 



Hq = Etc -h^^(Tf -K 
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where b runs over the bonds of a square lattice and Htc is 
given by Eq. ([T]). Note that the fields and induce 
different types of phase transition. The term with hz 
creates virtual pairs of Z2 charges, which condense when 
the field strength exceeds a certain threshold. This phe- 
nomenon may be described as a Higgs transition, or as 
vortex confinement. By duality, the field causes the 
condensation of vortices and charge confinement. The 
competition of the two terms should result in an inter- 
esting phase diagram, which is the subject of this paper. 

We approach the problem by reducing the quantum 
Hamiltonian to a classical anisotropic Z2 gauge Higgs 
Hamiltonian on a three-dimensional cubic lattice; see 
Eq. ^ below. We expect the phase diagram to be qual- 
itatively similar to that for the isotropic case, i.e., model 
M^,2 as defined by Wegner [§]. Monte-Carlo simulations 
have been performed for the latter model because it is 
more amenable to numerics. Some properties of the 
phase diagram in the isotropic case were predicted by 
Fradkin and Shenker In particular, the topological 
phase is bounded by second-order lines corresponding to 
charge condensation (for hz Jx^Jz) and vor- 

tex condensation (for hz ^ hx ~ Jx,Jz), but the two 
condensate phases are continuously connected. For the 
quantum Hamiltonian ([2]), a connecting path is realized 
by increasing hz so as to polarize the spins in the z di- 
rection, rotating the field in the xz-plane, and decreasing 
it again. However, the two phase transitions are clearly 
different, therefore the corresponding lines cannot join 
smoothly. 

A previous numerical study involving 10'^ sites by 
Jongeward, Stack, and Jayaprakash [lo| showed the two 
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lines merge into a first-order line that partially separates 
the charge and vortex condensates, and suggested the 
2nd-order lines might become Ist-order before merging. 
As discussed below, our data for larger systems (up to 
GO'^ sites) do not confirm this conjecture. Thus, the 
point where all three lines meet is likely to be a novel 
type of multicritical point. Note that each of the 2nd- 
order transitions can be characterized by an Ising-type 
order parameter, i.e., the amplitude of the corresponding 
condensate. The two parameters must somehow coexist 
though they are incompatible at the classical level due 
to the nontrivial braiding between charges and vortices, 
(ie., the kinetic terms describing the charge and vortex 
transport do not commute.) 

The reduction to a classical problem. The Hamilto- 
nian Q is not gauge- invariant, but it can be mapped to a 
Z2 gauge theory by introducing a dummy spin variable /Xs 
( "matter field" ) at every site, but only considering states 
l^") such that Hsl'^) = l^*). This constraint is a proto- 
type of the gauge- invariance condition, /if Agj^) = l^*), 
which says that flipping the spin on a site and all incident 
bonds docs not change the state. To turn one constraint 
into the other, we apply the transformation: 
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where Uuv belongs to the bond connecting sites u and v. 
Then, the Hamiltonian becomes: 

7? = - J, ^ - J, ^ i?p - /i, ^ 

s p b 

-hz^ti-Kytil- (4) 

uv 

Note that in the first term we have replaced Ag by /is 
using the gauge- invariance condition, /tf^s = 1- 

We now map this 2-d quantum Hamiltonian onto a 
(2+l)-d classical one. The overall scheme is standard 
but some care should be taken to preserve the 
gauge invariance. We let At = f3/n, and approximate 
the quantum partition function Z ~ Tr[exp(— /3_ff)'P] by 
Tr[exp(-AriJ^)7'cxp(-AriJ2)]", where V is the pro- 
jector onto the gauge-invariant subspace and , arc 
the terms in the quantum Hamiltonian that depend on 
cr^ , /if and (Tj , /ig , respectively. This expression can be 
written as a classical partition function on a cubic lat- 
tice. The classical variables ah,tilJ^s,t = il, in each time 
slice i, correspond to cr^ and /i^ respectively. But when 
we change from 2d to 3d, new vertical bonds (along the 
time direction) appear. The classical spins on the vertical 
bonds between two time slices correspond to a choice of 
term in the expansion of 7-" = Hs + A*s^s))- Thus 
we arrive at this classical Hamiltonian: 

Hc^-Y^ Ai|;-^^/i„(7„„/i„ - XI '^pi"^ n '^J ! 
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^bond = -2lntanh 

'^bond = 



A , = — Intanh/i,, 
pi 2 
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— vertical bonds; (6a) 

— horizontal bonds; (6b) 

— vertical plaquettes; (6c) 

— horizontal plaquettes, (6d) 



where Jx = JxAt, Jz = JzAr, — h^Ar, ~ hzAr. 
This model is an anisotropic generalization of the Z2 
gauge Higgs model Q . 

As a final step, we eliminate the redundancy by fixing 
/is. This only changes the classical partition function by 
a constant factor since Hamiltonian ([5]) can be written in 
terms of the gauge-invariant variables Suv = fJ-uC^uvfJ-v- 
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More detailed calculations show that the quantum and 
classical partition functions are related by 



Z = (isinh(2J^; 



k/2 / ^ \ 7x1/2 

(isinh(2/i,)J Zc, (8) 



where k and m are the number of vertical bonds and 
plaquettes, respectively. 

Of course, Eq. ([8]) holds only in the limit Ar 0. 
However, we take the liberty of parametrizing the general 
classical Hamiltonian ^ by Jx, Jz, hx, hz, even though 
the corresponding quantum problem may not be defined. 
In the isotropic case, two parameters will suffice: 
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where Abond = hz, Xpi 
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ilntanh/ij:. This model is 



equivalent to the isotropic Z2 gauge Higgs model Q and 
in what follows we compute its phase diagram. 

Phase diagram in the isotropic case. At Abond ~ all 
configurations, including ground states, have the same 
degeneracy factor 2^^. The actual physical variables 
in this limit are plaquette numbers Np = Yij^p^jy ^^'^ 
the model itself is dual to the 3D classical Ising model 
(Eq. ^ is also known as the 3D Ising gauge theory [^). 
Using high-accuracy results of Ref. [l^l for the critical 
point and the duality relation Api = — 1/2 lntanh(J/r), 
where J is the Ising exchange coupling, we obtain Ap'j'' = 
0.7614125. 

At arbitrary values of Abond and Api the model is self- 
dual [ist . i.e. it maps to itself under the coupling con- 
stant transformation Abond, pi ^ ^ 1/2 Intanh(Api^bond)- 
This means that the phase diagram has a symmetry, or 
self-duality, line defined by Abond = — 1/2 Intanh(Api). 
Under the duality mapping (Abond = 0, Api = 
0.7614125) (Abond = 0.221655, Api = 00), which gives 
us two Ising-type critical points on the phase diagram. 



3 



To calculate the rest of the phase diagram we per- 
formed Monte Carlo simulations using standard single- 
spin flip updates, supplemented by rare (once per A''^ up- 
dates) flips of all spins belonging to bonds cut by planes 
oriented along any one of the crystal axes, or along any of 
the diagonals to these axes. There are 9N possible planes 
satisfying this condition, and we select any of them at 
random. The plaquette energy (second term in ^) is 
conserved by this update. To determine the 2nd-order 
critical lines, we employed a standard finite-size scaling 
analysis of the specific heat Cy, for linear system sizes 
N = 24, 36, 48, and 60 (ic., for spins). First-order 
critical points were identified and located using energy 
distributions. These distributions are bi-modal (have two 
maxima) for the first-order transitions and single-modal 
otherwise. We thermalized our samples for up to 10^ MC 
sweeps (one sweep having elementary updates). The 
data were accumulated for ^ 4 x 10^ MC sweeps. 

The resulting phase diagram is presented in Fig[T] The 
first-order transition coinciding with the self-duality line 
was observed for 0.2575(5) > Abond > 0.22635(5). Out- 
side of this interval we saw no bi-modal structure in the 
energy distribution for system sizes up to = 60. The 
inset of Figl2] shows the evolution of the energy distri- 
bution function along the self-dual line. Even when the 
bi-modal structure is observed it is extremely weak, de- 
veloping only for large N, and the distribution can be 
sampled in the minimum without flat-histogram or sim- 
ilar reweighting techniques. 
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FIG. 1: (color online). The phase diagram of the Hamilto- 
nian (|9]). Circles correspond to the second-order transitions 
(open and filled symbols are related by the duality transfor- 
mation). Filled squares describe the first-order self-dual tran- 
sitions. Bold and dashed lines are used to guide an eye and 
correspond to the 1st- and 2nd-order transitions, respectively. 
The phases are: (I) - topological phase; (II) - topologically dis- 
ordered phase; (III) - magnetically ordered phase. In inset (a) 
we show the region where all phases meet each other. In inset 
(b) we show three alternative ways of connecting the lines. 

As noted above, these results conflict with previous 



MC simulations in Ref. [lO[, who suggest the Ist-ordcr 
line splits into two Ist-order lines. The inset (a) of Fig[T] 
shows a closeup of the controversial region. Though we 
were able to resolve critical points with an accuracy of at 
least three digits, we observed no splitting of the self-dual 
Ist-order line into two Ist-order transitions. We also find 
no evidence for tri-critical points on the Ising-type lines 
as long as we can resolve two separate transitions. There 
remains a tiny parameter range between the apparent 
disappearance of the bi-modal distribution on the self- 
dual line (this disappearance probably due to our limited 
system size) and two resolved 2nd-order transitions. 
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FIG. 2; (color online). The distance between two maxima 
in bi-modal energy distributions along the self-duality line for 
TV = 60 as a function of Abond. The inset shows examples of 
the energy distributions at various values of Abond- 

To probe the behavior in this tiny parameter range we 
need a different approach. We therefore scanned energy 
distributions at 30 points (A^ = 48) along the line perpen- 
dicular to the self-duality line right in the questionable re- 
gion (short solid line in the inset (a)). If the Ist-order line 
were to split above the scan, the third maximum would 
have to emerge in the energy distribution right between 
the two maxima we observe on the self-dual line - imply- 
ing that the energy maxima on the self-dual line could not 
merge smoothly, and right below the split, three maxima 
would have to be seen in the energy distribution. How- 
ever all distributions along the scan were found to have 
only one peak. It is also clear from the main part of 
Figl2] that on the self-duality line, the energy maxima 
approach each other and merge continuously as Abond in- 
creases. The curves presented in Figl5] follow a power 
law near the vanishing point, with corresponding critical 
exponent ~ 0.55. 

We thus conclude that the split Ist-order scenario does 
not work. Instead there arc three possibilities. Either 
all three lines merge at one point (case (1) in the in- 
set (b), Fig[T|); or the Ist-order line ends before or after 
the point where two 2nd-order lines touch the self-dual 
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line (cases (2) and (3) in the inset (b), FigH]). Unfor- 
tunately our data cannot distinguish between the alter- 
natives because the 2nd-order lines seem to touch at ex- 
tremely small (possibly zero) angle. Formally, option (2) 
fits the data best. Theoretically, the last two scenarios 
are less demanding since they fit the existing theory of 
phase transitions (our data suggest that the 2nd-order 
transitions cannot merge into a single smooth curve and 
form a kink at the self-dual line). We arc not aware of 
any effective theory leading to the first scenario. 
Phases. Using the two correspondence equations 

/iz = -ilntanh(ii,) = ^bond; (10a) 
Jz = -ilntanh(/i^) = Api, (10b) 

we can reformulate the phase diagram FigH] in terms of 
the renormalizcd parameters J^-, Jz, and hz of the 
TCM. The resulting phase diagram in terms of the ex- 
ternal fields is presented in FiglJl Let us go through the 
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FIG. 3: (color online). The phase diagram Fig[l]in terms of 
the renormalized external fields, hx = hx/S.T and hz = hz/S.T, 
of the model ((S} . The phases (I) , (II) and (II) are the same 
as in Fig[T] 

The phase (I) corresponds to the topological phase of 
the model (the "free charge" phase of the isotropic 
Z2 Higgs model (HM), [^). In this phase the system 
tends to have all Bp = 1 and = \ and a realization 
of such a state is obviously not unique. The plaque- 
ttcs with Bp = — 1 (magnetic vortices) and vertices with 
As = —1 (electric charges) appear mainly in the vicin- 
ity of the critical lines between the phases (I) and (II) 
(vortices) and (I) and (III) (charges). The phase (III) 
may be called "magnetically ordered" since the spins are 
mostly polarized in the z-direction. However, (cr^) also 
has nonzero value everywhere in the phase diagram. The 



true order parameter may be written as (/z^) using the 
gauge-symmetrized Hamiltonian (|4]). A non-zero value 
of this parameter results in the confinement of magnetic 
vortices (no free vortices) and the condensation of elec- 
tric charges. In the HM this is the "Higgs" phase. The 
phase (II) is characterized by a dual order parameter 
related to (u^), which can be defined by rewriting the 
Hamiltonian in different variables. Its nonzero value re- 
sults in non-conservation of total magnetic "charge" and 
condensation of magnetic vortices while electric charges 
are confined (no free charges). This phase corresponds 
to the "confinement" phase of the HM. The transition 
between the phases (II) and (HI) is accompanied by a 
sharp change in the number of vortices and charges, cor- 
responding to a "liquid-gas" type transition. The self- 
duality symmetry reflects the symmetry between charges 
and vortices. 

Summary. The topological phase of the toric code 
model (the "free charge" phase of the 3d gauge Higgs 
model) remains stable in a rather wide range of fields 
and breaks down via two Ising type transitions whose 
critical lines meet with the Ist-ordcr one corresponding 
to a liquid-gas type transition. The Ist-order line either 
meets with two 2nd-order lines in one multicritical point, 
or terminates before or after the point where two 2nd- 
order lines touch the self-duality line. The construction 
of an effective field theory for this multicritical region is 
an interesting open problem. 

We thank E. Fradkin, B. Svistunov, S. Trebst, M. 
Troyer, I. Affleck, and K. Shtengel for discussions. We 
are also indebted to M. Berciu and J. Heyl whose research 
clusters were used to perform our MC simulations. 



[1] P.W.Anderson, Science 235, 1196, (1987). 

[2] N. Read, B. Chakraborty, Phys. Rev. B 40, 7133 (1989). 

[3] A.Yu. Kitaev, Ann. Phys. (N.Y.) 303, 2 (2003). 

[4] G. Misguich, D. Serban, V. Pasquier, Phys. Rev. Lett. 

89, 137202 (2002); cond-mat/0204428. 
[5] N. Read, S. Sachdev, Phys. Rev. Lett. 66, 1773 (1991). 
[6] R. Moessner, S. L. Sondhi, Phys. Rev. Lett. 86, 1881 

(2001); cond-mat/0007378. 
[7] S. Trebst, P. Werner, M. Troyer, K. Shtengel, and Ch. 

Nayak, PRL 98, 070602 (2007). 
[8] F.J. Wegner, J. of Math. Phys. 12, 2259 (1971). 
[9] E. Fradkin and S. Shenker, Phys. Rev. D 19, 3682 (1979). 
[10] G.A. Jongeward, J.D. Stack and C. Jayaprakash, Phys. 

Rev. D 21, 3360 (1980). 
[11] M. Suzuki, Progress on Theor. Phys. 56, 1454 (1976). 
[12] R. Gupta and P. Tamayo, The critical exponents for the 
3d Ising model, US-Japan Bilateral Seminar - Maui, Au- 
gust 28-31, 1996. 
[13] R. Balian, J.M. Drouffe, and C. Itzykson, Phys. Rev. D 
11, 2098 (1975). 



